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1  Objectives 

1.  Time-frequency  filtering  of  jammed  GPS  signals; 

2.  Carrier-phase  ambiguity  resolution  in  the  presence  of  noise. 

2  Description  of  Achievements  of  Objectives 

Both  Objectives  have  been  accomplished.  Results  have  been  published  in  two  refereed-journal  papers  and  two 
refereed  conference-proceeding  papers.  Two  final  papers  are  being  prepared  for  submission. 

2.1  Time-frequency  filtering  of  jammed  GPS  signals 

For  Test  and  Evaluation  (T&E)  systems  on  DOD  test  ranges,  GPS  has  been  used  as  the  primary  source  of  Time- 
Space-Position  Information  (TSPI).  The  accuracy  requirement  for  TSPI  is  in  the  sub-meter  range,  an  order  of 
magnitude  higher  than  that  required  for  actual  combat  systems.  To  achieve  the  high  accuracy,  it  is  necessary  to 
use  carrier-phase  tracking  of  GPS  signals.  To  emulate  the  characteristics  of  enemy  GPS  threats  in  the  battlefield, 
the  TSPI  test  environment  also  includes  GPS  jamming  sources. 

The  spread-spectrum  structure  of  GPS  signals  provides  inherent  jamming  tolerance  for  GPS  receivers.  After 
despreading  the  received  signal,  the  original  satellite  signals  are  collapsed  into  a  narrow  band  around  the  carrier 
frequency,  while  the  jamming  signals  are  spread  due  to  the  lack  of  correlation  with  the  PN  codes  that  encode 
the  satellite  information.  The  portion  of  the  spread  jamming  signal  remaining  within  the  frequency  band  of  the 
tracking  loop  enters  the  loop  together  with  the  satellite  signal.  Thus  signal  tracking  can  be  stable  only  if  the 
jamming-to-signal  ratio  (JSR)  is  below  the  processing  gain  of  the  spreading  code.  In  general,  a  JSR  of  greater 
than  40dB  will  prevent  the  GPS  receiver  from  tracking  the  satellite  signals  and  estimating  its  own  position.  In  a 
hostile  environment  where  jamming  sites  may  be  close  to  GPS  users,  a  larger  JSR  is  possible.  How  to  achieve  the 
desired  accuracy  for  a  GPS-based  TSPI  system  in  the  presence  of  strong  jamming  is  an  important  but  outstanding 
problem. 

In  the  project  period,  we  developed  a  class  of  time-frequency  filters  based  on  the  combination  of  the  empirical¬ 
mode  decomposition  (EMD)  method  and  a  general  blind-source  separation  (BSS)  algorithm.  We  obtained  evi¬ 
dence  that  the  method  is  able  to  separate  jamming  from  the  GPS  signal  for  JSR  up  to  45dB. 

A  forefront  research  area  in  signal  processing  is  particle  filters.  The  idea  is  to  evolve  the  probability  distri¬ 
bution  of  a  signal  by  using  a  large  number  of  particles  according  to  the  system  equations  and  some  stochastic 
processes,  which  is  similar  to  Monte-Carlo  simulation  in  physics  and  chemistry.  Motivated  by  the  fact  that  par¬ 
ticle  filters  have  been  used  widely  in  various  types  of  signal-processing  tasks,  we  applied  this  technique  to  GPS 
positioning  of  moving  objects  in  a  jamming  environment.  In  particular,  we  considered  a  class  of  regularized 
particle  filters,  suitable  for  estimating  the  position  of  a  moving  object  (e.g.,  a  car)  equipped  with  some  proper 
GPS  C/A  code  receiver.  Theoretically,  a  question  of  interest  is  how  the  estimation  error  depends  on  uncertainties 
in  the  velocity  measurement  of  the  car  and  on  the  noise  level  in  the  GPS  signal.  Our  analysis  of  the  covariance 
matrix  constructed  from  simulated  particles  led  to  a  formula  relating  this  matrix  to  the  covariance  matrices  of  the 
velocity  and  of  the  position  error  from  least-squares  processing  of  GPS  pseudoranges.  The  formula  was  verified 
by  numerical  simulations. 

Practically,  in  order  to  address  the  problem  of  occasional  but  inevitable  large  errors  (outliers)  in  the  GPS 
observations,  we  developed  a  robust  particle-filtering  technique.  We  demonstrated  that  the  strategy  is  effective 
for  mitigating  the  effect  of  outliers  for  both  Gaussian  and  non-Gaussian  noise  sources.  Even  in  the  absence  of 
outliers,  our  strategy  can  be  useful  for  improving  the  positioning  accuracy. 
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•  V.  Kamath,  Y.-C.  Lai,  S.  Urval,  and  L.  Zhu,  “Empirical  mode  decomposition  and  blind  source  separa¬ 
tion  methods  for  antijamming  with  GPS  signals,”  pp.  335-341  in  IEEE  PLANS  (Position  Location  and 
Navigation  Symposium),  April  24-26,  2006,  San  Diego,  CA. 

i 

•  L.  Huang  and  Y.-C.  Lai,  “Sequential  Monte  Carlo  scheme  for  Bayesian  estimation  in  the  presence  of  data 
outliers,”  Physical  Review  E  75,  056705(1-6)  (2007). 

•  V.  Kamath  and  Y.-C.  Lai,  “Empirical  mode  decomposition  with  applications  to  noise  reduction  and  anti¬ 
jamming,”  in  preparation. 

2.2  Carrier-phase  ambiguity  resolution  in  the  presence  of  noise 

In  order  to  determine  a  GPS  receiver’s  position,  the  following  two  pieces  of  information  are  needed:  (1)  the 
satellite  positions  and  (2)  the  traveling  times  of  signals  or  the  numbers  of  wavelength.  With  these,  the  receiver’s 
position  can  be  calculated  by  a  simple  trilateration  procedure. 

Information  about  the  satellite  positions  is  encoded  in  the  GPS  signal  via  the  traditional  spread-spectrum 
method.  The  respective  distances  between  satellites  and  the  GPS  receiver  can  be  obtained  by  figuring  out  the 
numbers  of  wavelengths  of  the  carrier  signals.  The  fractional  parts  of  the  wavelength  can  be  measured  by 
phase-locked  loop  based  detectors  embedded  in  the  receiver.  The  integer  parts  are  ambiguous  and  these  inte¬ 
ger  ambiguities  must  be  resolved  accurately  to  meet  the  T&E  requirement  of  sub-meter  precision  in  a  jamming 
environment.  There  exist  several  integer-parameter  estimation  methods,  but  how  they  perform  under  strong,  non- 
Gaussian  noise  or  jamming  is  unknown.  During  the  project  period,  we  evaluated  the  performance  of  a  state-of-art 
integer-parameter  estimation  algorithm  and  explored  improvements  for  use  in  a  jamming  environment. 

Particle  filters  are  naturally  suitable  for  integer-parameter  estimation  with  GPS  carrier  waves,  as  the  system 
equations  can  be  obtained  via  straightforward  geometrical  arguments.  During  the  project  period  we  successfully 
demonstrated  the  feasibility  of  this  approach. 

Publications 

•  M.  Shah  and  Y.-C.  Lai,  “Performance  of  integer  parameter  estimation  algorithm  for  GPS  signals  in  noisy 
environment,”  pp.  166-174  in  ION  GNSS  1 7th  International  Technical  Meeting  of  the  Satellite  Division , 
Sept.  21-24,  2004,  Long  Beach,  CA. 

•  L.  Zhu,  Y.-C.  Lai,  M.  Shah,  and  S.  Mahmood  “Efficiency  of  carrier-phase  integer  ambiguity  resolution  for 
precise  GPS  positioning  in  noisy  environments,”  Journal  of  Geodesy  81,  149-156  (2007). 

•  A.  Rammohan,  L.  Huang,  and  Y.-C.  Lai,  “Integer  ambiguity  resolution  using  particle  filters,”  in  prepara¬ 
tion. 


3  New  Findings 

3.1  Performance  of  integer  least-squares  method  for  ambiguity  resolution  with  GPS  signals  in 
the  presence  of  noise 

3.1.1  Background 

In  general,  centimeter  GPS  positioning  accuracy  requires  a  precise  tracking  of  the  carrier  phase  that  consists 
of  two  parts:  a  directly  measured  fractional  part  (with  measurement  error  at  millimeter  level)  and  an  unknown 
integer  part,  also  called  the  integer  ambiguity.  The  key  to  carrier-phase-based  precise  positioning  is  to  resolve  the 
integer  ambiguity,  which  is  an  extremely  challenging  task,  particularly  when  large  noise  or  jamming  is  present. 
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Existing  ambiguity  resolution  techniques  can  be  divided  into  several  categories  [1].  The  first  category  in¬ 
cludes  the  simplest  techniques  which  use  C/A-code  or  P-code  pseudoranges  directly  to  determine  the  ambiguities 
of  the  corresponding  carrier  phase  observations.  The  precision  of  the  code  range  is  not  good  enough  to  deter¬ 
mine  the  integer  ambiguities  and  inter  frequency  linear  combinations  are  usually  used  for  estimating  ambiguities 
[2].  The  second  category  of  algorithms  employ  the  Ambiguity  Function  Method  (AFM)  [3].  This  technique 
uses  only  the  fractional  value  of  the  instantaneous  carrier-phase  measurement  and,  hence,  the  ambiguity  func¬ 
tion  values  are  not  affected  by  the  whole  cycle  change  of  the  carrier  phase  or  by  cycle  slips  (also  see  [2,  4]). 
The  third  category  comprises  the  most  abundant  group  of  techniques  which  are  based  on  the  theory  of  integer 
least  squares  [5,  6,  7,  8,  9,  10,  11,  12,  13,  14,  15,  16,  17,  18].  Parameter  estimation  in  theory  is  carried  out  in 
three  steps:  float  solution,  integer-ambiguity  estimation,  and  fixed  solution.  Each  technique  makes  use  of  the 
variance-covariance  matrix  obtained  at  the  float  solution  step  and  employs  different  ambiguity  search  processes 
at  the  integer  ambiguity  estimation  step.  Based  on  certain  search  criterion  [19],  the  search  algorithm  can  utilize 
the  traditional  techniques  of  mathematical  programming  to  guide  the  global  optimization  [20,  12,  13]  and/or 
decorrelation  techniques  to  reduce  the  search  space  [9,  16,  14].  Guided  random  searching  techniques  can  be  used 
to  combat  nonlinearity  [21,  17].  Note,  however,  that  decorrelation  would  help  speed  up  searching  for  the  integer 
solution  only  if  the  dimension  is  not  too  large  [14]. 

For  any  method,  a  practical  concern  (especially  for  kinematic  positioning)  is  that  the  resolution  time  of  the 
integer  ambiguity  is  sensitive  to  the  carrier-phase  measurement  noise.  In  a  noisy  environment,  e.g.,  battle  field 
with  strong  jamming  signals,  radio  frequency  interference  (RFI)  signals  are  spread  in  the  frequency  domain  by 
the  de-spreading  process.  These  spectrally  spread  RFIs  increase  the  effective  noise  floor  in  a  GPS  receiver,  mak¬ 
ing  carrier-phase  measurements  more  noisy  and  hence  the  time  to  resolve  integer  ambiguity  longer.  Whenever  a 
loss  of  track  due  to  jamming  or  blockade  occurs,  the  integer  ambiguity  has  to  be  resolved  again,  during  which 
no  high-precision  position  information  can  be  available.  To  reduce  the  phase-measurement  noise  and  the  integer- 
ambiguity  resolution  time,  various  techniques  such  as  those  based  on  signal  processing  (e.g.,  time-frequency 
domain  processing)  [22],  subspace  processing  [23]  and/or  receiver  antenna  design  (e.g.,  beam  forming,  null 
steering)  have  been  proposed  [24,  25]  for  antijamming.  However,  even  with  the  application  of  these  techniques, 
typically  there  are  still  residual  errors  in  the  phase  measurement  and  greater  efforts  are  required  for  fast  reso¬ 
lution  of  the  integer  ambiguity.  Given  a  minimal  requirement  for  positioning  accuracy  and  acquisition  time,  it 
is  desirable  to  know  the  corresponding  noise  level  of  carrier-phase  measurements,  in  order  to  employ  only  the 
necessary  jamming  mitigation  techniques. 

During  the  project  period,  we  investigated  the  performance  of  a  typical  integer  least-squares  ambiguity  reso¬ 
lution  algorithm  in  noisy  environments.  In  particular,  we  addressed  how  the  convergence  (acquisition)  time  of  the 
algorithm  depends  on  the  carrier-phase  measurement  noise.  Naively,  one  may  expect  that  the  time  and  the  noise 
variance  a2  have  a  linear  relationship,  as  suggested  by  the  fact  that  the  variance  of  n  independent  noisy  mea¬ 
surements  of  variance  a2  is  a 2 In.  However,  our  theoretical  analysis  showed  that  the  acquisition  time  depends 
linearly  on  the  standard  deviation  of  noise  (or  the  noise  amplitude)  cr,  not  on  a2,  as  also  verified  by  numerical 
simulations  using  a  generic  integer-parameter  estimation  algorithm  [16].  Since  the  naive  relationship  n  oc  a2 
would  require  much  more  observation  samples  than  the  linear  relation  n  oc  a  (for  large  a  and  n),  the  finding  is 
important  for  the  design  of  precise  GPS  positioning  system  in  noisy  environments  with  tight  constraint  on  time, 
such  as  in  kinematic  GPS  positioning.  Our  finding  is  particularly  encouraging  as  it  suggests  the  possibility  of 
having  integer-parameter  estimator  to  achieve  significantly  short  convergence  time  for  precise  positioning. 

In  the  following,  we  briefly  review  a  standard  GPS  model  and  describe  a  generic  ambiguity  resolution  algo¬ 
rithm  based  on  integer  least  squares.  Our  focus  is  on  the  performance  of  the  algorithm  in  the  presence  of  noise. 
We  shall  also  derive  a  theoretical  relationship  between  the  convergence  time  and  the  noise  amplitude  and  provide 
numerical  verification. 
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3.1.2  System  model 


GPS  signals  are  usually  corrupted  by  jamming  and  several  other  forms  of  errors.  These  error  sources  can  be 
partially  cancelled  by  using  the  technique  of  double-differencing.  For  positioning  based  on  the  carrier  phase,  we 
assume  that  jamming  and  the  residual  errors  can  be  modeled  as  noise  in  the  phase.  Under  this  consideration,  the 
distance  between  a  satellite  i  and  a  receiver  A  can  be  modeled  as 

p\(k)  =  \[$'A{k)  +  n'A]  +  v(k),  (1) 


where  k  is  the  discrete  sampling  time,  A  is  the  wavelength  of  the  carrier,  $  is  the  fractional  part  of  the  carrier 
phase,  n  is  the  integer  part  of  the  initial  carrier  phase  (i.e.,  at  k  =  0),  and  v  is  the  modeling  error.  The  range  p  can 
be  expressed  in  terms  of  the  satellite  position  [ xt(k)1yl(k),zt(k )]  and  receiver  position  [x^fc),  yA(k),  zA(k)]: 

p\{k)  =  {(x‘(fc)  -  xA(k )]2  +  [y'(k)  -  yA(k)]2  +  [z'(fc)  -  zA{k)]2}1/2.  (2) 


The  goal  is  to  calculate  the  receiver  position  [xA(k),  yA(k),  zA(k)]  by  using  known  positions  of  at  least  three 
satellites  and  the  corresponding  phase  measurements  <1>. 

To  linearize  Eq.  (2),  we  let  xA  =  xA0  +  Axa,  yA  =  yA0  4-  A yA,  zA  =  zA0  +  A zA,  where  [x^o,  yA0,  zA0] 
is  a  known  reference  position  near  the  receiver,  which  can  be  estimated  using  code  pseudoranging.  Substituting 
the  linearized  version  of  Eq.  (2)  into  Eq.  (1)  yields 


A(k )  -  pA0(k)  =  alx{k) Ax  +  aly(k)Ay  +  a'z(k)Az  -  \n\ , 


(3) 


where 


4(fc)  =  - 


X*(/c)  -  X/\Q 

Pa o(^) 


y'(k)  -  yM 

P\ o(k) 


and  a\{k)  —  - 


z'(fc)  -  2,40 

P\  o(fc) 


(4) 


Suppose  the  receiver  is  static  and  three  satellites  are  continuously  tracked  from  epoch  0  to  epoch  n,  the  linearized 
observation  equations  can  be  expressed  in  matrix  form  as 


y  =  Ax  +  Bz  +  v, 


(5) 


where 

y  =  [A<4(0)  -  4o(0),  (0)  -  p\ o(0), . .  ,]T, 


X  = 

by 


[Ax,  Ay,  A z]T,  z  =  [n\,  n2A,  n3A]Ty  v  is  the  measurement  noise  vector,  and  the  matrices  A  and  B  are  given 


A  = 


4(0) 

4(0) 

4(o)  1 

'  -A 

0 

0 

4(o) 

4(°) 

4(o) 

0 

-A 

0 

4(0) 

4(o) 

4(o) 

0 

0 

-A 

4(i) 

4(i) 

4(i) 

,  and  B  = 

-A 

0 

0 

4(D 

4(i) 

4(i) 

0 

-A 

0 

4(i) 

4(D 

4(i) 

0 

0 

-A 

: 

3(n+l)x3 

•  - 

3(n+l)x3 


(6) 


3.1.3  Integer  Least  Squares 

In  the  linearized  observation  equation  Eq.  (5),  x  and  z  are  unknown  real  and  integer  vectors,  respectively,  which 
are  to  be  estimated  from  the  satellite  data  by  using  the  maximum  likelihood  (ML)  estimator: 


(xa/l,za/l)  =  argmaxPy)xz(y|x,z), 

X  ,z  ’ 


(7) 
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where  (x,z)  E  Rp  x  Zq ,  and  py|x,  Z(y|x,z)  is  the  probability  of  observing  y  given  x  and  z.  Assuming  the 
stochastic  process  v  is  Gaussian  with  zero  mean  and  covariance  E,  we  have 

(xml,zml)  =  arg min(y  -  Ax  -  Bz)TE_1(y  -  Ax  -  Bz).  (8) 

x,z 

If  we  consider  a  block  matrix  [A  B]  and  the  unknown  vector  [xz]T,  the  floating  solution  [x  z]7  can  be  obtained 
by  solving 

[AB]tE-1[AB][xz]t  =  [A  B]rE_1y>  (9) 

Eliminating  x  leads  to  the  least-squares  solution  of  z: 

z  =  rBr(I  -  A(ATE_1  A)-1  Ar)y,  (10) 

where  T  =  (BrC'B)_1  is  the  covariance  matrix  for  z,  and  C'  =  E_1  -  E-1  A  (A7  E"1  A)-1  A7  E_1. 

In  general,  z  is  a  real  vector  and  the  ML  solution  can  be  found  in  the  following  least-squares  sense  [7,  12]: 

zML  =  argmin(z  -  z)Tr_1(z  -  z)  for  z  E  Zq .  (11) 

z 

If  r  is  diagonal,  the  solution  for  z  can  be  found  by  rounding  each  component  of  z  to  its  nearest  integer.  This 
simple  approach,  however,  does  not  work  in  realistic  situations  where  the  matrix  T  is  typically  not  diagonal.  As  a 
result,  simple  rounding  off  will  not  give  the  correct  estimate  of  z  and  the  search  for  the  true  integer  vector  has  to 
be  performed  over  the  entire  integer  space  by  using  some  efficient  searching  algorithm  [5,  6,  1 1 , 8,  1 6,  1 8].  Once 
zml  is  found,  the  estimate  of  the  real  position  vector  x  can  be  obtained  by  substituting  z ml  into  the  least-squares 
equation  (8).  This  yields 

xml|z  =  (ArE-1A)-1ArE-1(y  -  BzML).  (12) 

It  should  be  noted  that  the  ambiguity  fixing  process  (Eq.  1 1)  breaks  down  Eq.  (9),  resulting  in  the  above  least 
squares  ambiguity  search  (LSAS)  algorithm  not  being  an  optimal  one,  i.e.,  it  is  not  guaranteed  that  (xml\z^ml)  = 
(xml,zml)-  It  has  been  shown  that  ( *ml\z » z  ml)  can  be  away  from  (xml,zml)  and  a  more  general  searching 
criterion  has  to  be  used  [19,  4,  26]. 

t 

3.1.4  Scaling  relation  between  convergence  time  and  noise  amplitude 

In  practice,  we  are  most  interested  in  cases  when  the  positioning  accuracy  is  high,  or  the  probability  of  cor¬ 
rect  estimate  of  the  integer  ambiguity,  Pc ,  is  close  to  one.  Particularly,  we  wish  to  know,  given  certain  noisy 
environment,  how  long  it  will  take  us  to  achieve  a  desired  value  of  Pc,  where  Pc  is  close  to  one.  This  can 
be  addressed  by  investigating  the  influence  of  noise  on  Pc.  In  this  study,  Pc  obtained  from  LSAS  is  used  and 
(xa*l|z,za/l)  =  {xml,zml)  is  assumed  for  simplicity,  since  we  are  only  concerned  with  large  Pc  and  LSAS 
algorithm  generally  produces  good  results  in  this  situation. 

Assume  y  is  a  Gaussian  variable,  then  z,  a  linear  function  of  y,  is  Gaussian  too.  We  write  z  =  z  -f  u,  where 
u  is  Gaussian  with  zero  mean  and  covariance  matrix  I\  Multiplying  both  sides  by  G  =  T-1/2  and  defining 
y  =  Gz,  we  get  y  =  Gz  -f  u,  where  u  is  a  Gaussian  random  variable  with  zero  mean  and  unit  variance. 
Equation  (11)  can  then  be  written  in  an  equivalent  form 

2ml  =  argmin||y  -  Gz||2  for  z  €  Zq.  (13) 

The  set  {Gz|z  E  Zqj  constitutes  a  lattice  in  Rq.  Equation  (13)  suggests  that  the  maximum  likelihood  value  of  z 
can  be  found  by  computing  the  nearest  lattice  point  to  vector  y .  The  probability  Pc  that  zA/^  is  true  is  completely 
determined  by  the  Venoroi  cell.  By  definition  of  the  lattice,  a  lower  packing  distance  between  two  neighboring 
points  can  be  computed  by  using 

doc|G|1/<?,  (14) 
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where  |G|  is  the  determinant  of  the  matrix  G  and  q ,  the  number  of  unknown  integers,  is  the  dimension  of  G.  Asa 
rough  estimate,  the  distance  d  can  be  used  to  compute  the  lower  bound  of  the  probability  Pc.  Alternatively,  given 
a  desired  value  of  Pc,  there  is  a  corresponding  value  of  d  that  depends  on  factors  such  as  the  noise  amplitude  and 
the  number  of  data  samples.  We  emphasize,  however,  that  such  use  of  the  rough  estimate  of  the  lower  bound  of 
Pc  is  only  for  the  purpose  of  obtaining  a  scaling  dependence  of  the  convergence  time  on  the  noise  strength,  which 
by  its  nature  is  not  precise.  In  fact,  precise  estimates  of  the  lower  and  upper  bounds  have  been  obtained  recently 
[18]. 

To  be  able  to  find  a  closed-form  solution  for  |G|,  here  we  consider  a  two-dimensional  (2D)  GPS  setup  as 
shown  in  Fig.  1.  All  results  based  on  this  2D  model  can  be  extended  to  the  real  three-dimensional  GPS  setup. 
Assume  that  there  are  q  satellites  tracked  and  the  noise  terms  for  the  corresponding  carrier-phase  measurements 
are  independent  of  each  other  with  covariance  matrix  E  =  a2 1.  The  inverse  of  T  becomes 


where 


r_1  =  -^[BrB  -  BtA(A7'A)-1AtB]  =  ^ 

Gl  Gl 


nl  - 


d\d$  —  (& 


2  J 


n  q 


n  q 


n  q 


d\  =  Elax(fc)l2'  d2  =  ^^4(fc)a^(A:),  and  d3  =  ]T  ^[^(A:)]2, 


fc= 0  1=1 


k= 0  i=l 


k= 0 i= 1 


and  Q  is  a  q  x  q  matrix  with  elements  given  by 


(15) 


(16) 


n  n  n  n  n  n 

Qij  =  +  di  (17) 

k= 0  k= 0  fc=0  k= 0  A:=0  k=  0 

Apparently,  IGI1/9  =  |r— 1 1 !/2^  depends  on  both  g  and  n. 

To  find  the  relationship  between  g  and  n  for  a  given  PCy  or  equivalently,  a  constant  |r~ 1 1,  we  assume 
[xAO,yAo]  =  [0,0],  the  initial  angle  of  satellite  i  is  a,  the  orbital  speed  of  the  satellite  is  a,  and  the  sampling 
period  is  T.  We  then  have 


^  alx(k)  =  —  ^2  cos  (a  +  6tkT)  ~  ^  cos  a  +  sin  a  ^  akT  = 


k= 0 


T a  sin  a  —  (n  -h  1 )  cos  a, 


(18) 


7 


k 


where  higher  order  terms  of  akT  in  the  Taylor  expansion  are  neglected.  The  approximation  error  is  generally 
small  since  akT  is  relatively  small.  Similarly,  we  have 


Jt=0 

i>xWi2 

k= 0 

Ekmi2 

k= 0 
n 

Y^alx(k)aly(k) 
k= 0 


-n2Ta  cos  a  —  (n  - f  l)sina  (19) 

\n(n  -f  1)(2 n  +  l)T2a2  sin2  a  -  2n2Ta  sin  a  cos  a  4-  (n  4  l)cos2a  (20) 
6 

\n(n  4  1)(2 n  -I-  l)r2d2cos2  a  -  2n2Tasinacosa  4  (n  4  l)sin2  a  (21) 
o 

-n(n  4  1)(2 n  4  1)T2  a2  sin  a  cos  a  -  ]-n2Ta  4  (n  4  1)  sin  a  cos  a.  (22) 


Substituting  Eqs.  (18)  and  (19)  into  Eq.  (15),  one  can  see  that  each  entry  in  T_1  can  be  expressed  as  a  second- 
order  polynomial  of  n  divided  by  a2.  Thus,  for  increased  cr,  n  has  to  be  increased  proportionally  in  order  to 
maintain  a  desired  positioning  accuracy  (or  a  constant  |r_1|).  This  implies 


n  oc  a,  (23) 

for  n  »  1.  While  larger  noise  requires  more  data  samples  from  the  satellites,  the  linear  relation  indicates  that, 
to  achieve  a  desired  positioning  accuracy,  the  requirement  is  not  as  stringent  as  one  would  naively  expect  from 
n  oc  a2.  Our  derivation  suggests  that  if  the  satellites  were  kept  static  with  respect  to  the  static  GPS  receiver  [i.e., 
d  =  0  in  Eq.  (18)],  then  the  scaling  relation  n  oc  a2  would  hold.  It  is  the  relative  movement  between  GPS 
satellites  and  Earth  surface  that  introduces  another  degree  of  dependence  on  n,  resulting  in  a  shorter  estimation 
time  than  the  static  case.  Note  that  the  movement  of  the  receiver  will  have  similar  effect  if  the  velocity  of  the 
receiver  is  known  exactly.  This  is  in  fact  quite  encouraging  as  it  suggests  that  the  integer  parameter  estimation 
time  for  kinematic  GPS  is  likely  to  be  shorter. 


3.1.5  Simulation  Results 

Since  the  purpose  of  the  simulations  is  to  verify  the  theoretically  derived  dependence  of  the  covergence  time  on 
the  noise  amplitude,  it  is  necessary  to  be  able  to  vary  the  noise  variance  in  a  systematic  way.  A  synthetic  2D 
GPS  setup,  as  illustrated  in  Fig.  1,  is  suitable  for  this  purpose.  The  setup  is  similar  to  the  one  used  in  [16].  We 
assume  that  the  position  of  the  (GPS)  receiver  x,  which  is  to  be  determined,  can  be  modeled  as  a  zero-mean 
Gaussian  random  variable  with  certain  variance  in  each  dimension.  The  coordinate  axes  are  chosen  such  that 
the  origin  is  a  point  on  the  surface  of  the  earth  [a  point  on  the  periphery  of  a  circle  of  radius  equal  to  that  of  the 
earth  Re  =  6357km].  We  further  suppose  that  there  are  three  visible  satellites  orbiting  the  earth  at  the  altitude 
of  20,  200km  and  with  the  period  of  12  hours  (angular  velocity  of  l/120s_1).  The  satellites  transmit  a  carrier 
signal  of  wavelength  A  =  19cm  each,  and  their  coordinates  are  known  to  the  receiver.  The  receiver,  which  is 
assumed  to  be  completely  synchronized  with  the  satellites  (meaning  that  it  can  generate  the  transmitted  carrier 
signals),  measures  the  phases  of  the  received  carrier  signals  every  T  =  2s  and  unwraps  them  as  times  goes 
by.  By  multiplying  these  (unwrapped)  phase  measurements  by  the  wavelength  divided  by  2n  ,  the  receiver  can 
measure  its  distance  (or  range)  to  each  satellite  up  to  some  additive  noise,  which  is  assumed  to  be  yV(0.  a2)  and, 
of  course,  up  to  an  integer  multiple  of  the  wavelength.  (This  integer  multiple  can  be  thought  as  the  number  of 
carrier  signal  cycles  between  the  receiver  and  the  satellite  when  the  carrier  signal  is  initially  phase  locked.)  As 
we  have  described,  by  linearizing  the  range  equations,  the  problem  becomes  one  of  estimating  a  real  parameter 
x  (the  coordinates  of  the  receiver)  and  an  integer  parameter  z  (the  integer  multiples  of  the  wavelengths)  in  a 
linear  model.  In  the  simulation  that  follows,  the  actual  location  of  the  receiver  is  x  =  [50;  100]T,  which  will  be 
estimated  using  the  carrier-phase  measurements.  We  assume  that  the  standard  deviation  of  x  is  100m  along  each 
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Figure  2:  Pc  versus  time  for  a  =  1cm.  Note  that  Pc  after  t  =  200s  is  almost  equal  to  one  in  the  simulation. 


coordinate  axis.  The  satellites  make  angles  of  100,  130,  and  50  degrees  with  the  axis  initially,  and  the  direction 
of  rotation  for  all  of  them  is  clockwise.  The  standard  deviation  of  phase-measurement  noise  in  units  of  length  is 
assumed  to  be  in  the  centimeter  range.  Using  the  carrier-phase  measurements,  the  receiver  tries  to  find  its  own 
position  x  (as  well  as  the  ambiguous  integer  multiples  of  the  wavelengths)  as  a  function  of  time  by  solving  for 
the  ML  estimates. 

Figures  2  and  3  show  the  results  of  the  algorithm  for  a  =  1cm  in  terms  of  Pc  and  receiver  positioning 
respectively.  The  exact  value  of  the  Pc  is  computed  by  Monte  Carlo  simulation  of  1500  runs.  It  can  be  seen  that 
the  position  estimation  error  reduces  as  the  probability  of  correct  integer  estimation  approaches  unity.  When  the 
integer  ambiguity  is  resolved  correctly,  the  position  estimation  error  is  of  the  order  of  millimeters. 

Having  established  the  accuracy  of  the  integer  least-squares  algorithm  for  static  GPS  positioning,  we  wish 
to  evaluate  its  performance  in  the  presence  of  noise.  For  this  purpose,  the  algorithm  was  provided  with  various 
values  of  phase  measurement  noise  variance  as  input.  The  range  of  a  varied  from  0.2cm  to  2cm  with  increment 
of  0.2cm.  The  values  of  Pc  for  all  cases  were  calculated  using  1500  Monte  Carlo  simulations.  Figure  4  shows  Pc 
versus  time  for  various  o  values.  It  is  observed  that  as  a  increases,  the  time  it  takes  to  achieve  a  certain  level  of 
Pc  also  increases.  The  family  of  curves  in  Fig.  4  can  be  used  to  calculate  the  maximum  allowable  noise  variance 
for  a  given  amount  of  observation  time  and  required  value  of  Pc.  For  example  in  Fig.  4,  if  the  observation  time  is 
150  s  and  the  integer  ambiguity  estimate  is  to  be  reliable  with  90%  accuracy  (Pc  =  0.9),  the  maximum  allowable 
noise  variance  is  about  1.2cm. 

Figure  5  shows  the  linear  relation  characterizing  the  sensitivity  of  performance  to  noise  amplitude,  as  pre¬ 
dicted  by  our  theoretical  analysis.  A  use  of  this  result  is  that  the  performance  of  integer  least-squares  algorithm 
can  be  predicted  for  a  given  noise  level. 

3.1.6  Discussions 

We  have  addressed  the  performance  of  integer  least-squares  algorithm  for  GPS  signals  in  noisy  environments. 
Mathematically,  integer  ambiguity  resolution  is  equivalent  to  searching  for  the  integer  vector  closest  to  a  given 
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Figure  3:  Position  estimation  error  (meters)  versus  time  (sec),  where  the  initial  errors  (t  <  20s)  are  large  and  not 
shown  in  the  figure. 


Figure  4:  Pc  versus  time  for  a  values  ranging  from  0.2cm  to  2cm  (from  left  to  right).  Maximum  allowable  noise 
variance  for  given  values  of  time  and  Pc  can  be  determined  accordingly. 
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Figure  5:  Sensitivity  of  performance  to  measurement  noise  variance.  From  top  to  bottom:  Pc  = 
0.99, 0.95, 0.9, 0.75, 0.5  respectively. 


real  vector  on  an  integer  lattice.  In  a  noisy  environment,  the  probability  of  error  can  be  significantly  large.  We 
find  that  the  observation  time  required  to  achieve  a  fixed  value  of  a  lower  bound  of  Pc  and  thus  Pc  itself  is 
directly  proportional  to  the  standard  deviation  of  phase-measurement  noise,  in  contrast  to  the  naive  expectation 
that  the  time  is  proportional  to  the  variance  of  the  phase  noise.  This  suggests  the  possibility  of  achieving  short 
convergence  time  even  if  large  noise  is  present. 

It  should  be  noted  that  we  have  assumed  an  ideal  system  model  in  this  study.  For  example,  systematic  errors 
are  assumed  to  be  eliminated  by  using  double  differencing,  and  the  resulted  system  errors  are  approximately 
Gaussian  with  zero  mean.  This  may  not  be  true  in  practice,  especially  when  long  baselines  are  utilized.  In  this 
situation,  additional  modeling  of  large  residual  errors  has  to  be  employed,  e  g.,  the  means  of  residual  ionospheric 
and  tropospheric  errors  have  to  be  estimated  together  with  x  and  z.  If  the  means  of  these  residual  errors  remain 
constant  or  change  slowly  during  the  process  of  ambiguity  resolution,  the  convergence  time  of  ambiguity  reso¬ 
lution  will  only  be  delayed  approximately  by  a  constant,  the  overall  linear  relationship  between  the  time  and  the 
noise  amplitude  should  still  hold. 

We  also  remark  that  the  main  objective  of  our  work  is  to  obtain  a  scaling  relation  between  the  convergence 
time  and  the  noise  strength,  which  is  an  order-or-magnitude  type  of  estimate.  The  reason  is  that  the  convergence 
time  depends  on  too  many  algorithmic  and  system  details  and  it  is  not  possible  to  give  precise  numbers  to 
characterize  it.  For  this  purpose  we  have  based  our  theoretical  analysis  on  the  upper  and  lower  bounds  in  the 
probability  of  correct  resolution  of  the  integer  ambiguity  given  in  [16],  as  the  required  computational  time  is 
only  polynomial.  However,  these  bounds  can  be  poor  and  much  tighter  bounds  can  be  found  in  Ref.  [18].  We 
wish  to  emphasize  that  the  results  obtained  for  the  performance  in  the  presence  of  noise  are  based  on  the  integer 
least-squares  principle  in  general  and  thus  they  should  not  be  dependent  on  the  specific  search  method  used. 

For  future  directions,  parallel  algorithms  taking  advantage  of  multiple  satellite  signals  to  reduce  the  ambiguity- 
resolution  time  should  be  pursued.  In  addition,  the  study  for  the  effect  of  noise  on  GPS  positioning  should  be 
extended  to  kinematic  GPS  positioning  algorithms.  It  is  also  recommended  to  develop  a  particle-filter  based 
algorithm  and  to  compare  its  sensitivity  to  noise  with  those  of  integer  least-squares  algorithms  for  both  static  and 
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Figure  6:  Block  diagram  of  our  proposed  GPS  antijamming  scheme. 


kinematic  positioning. 

3.2  Empirical  Mode  Decomposition  and  Application  to  Anti-jamming  for  GPS  Signals 
3.2.1  Background 

GPS  signals  are  direct-sequence  (DS)  spread-spectrum  signals  and  therefore  have  some  degree  of  inherent  noise 
immunity  [27].  After  despreading  the  received  signal,  the  original  satellite  signal  can  be  collapsed  into  a  narrow 
bandwidth  about  the  carrier  frequency,  while  the  jamming  signals  are  spread  due  to  the  lack  of  correlation  with 
the  pseudo-random  code  that  encodes  the  satellite  information.  A  portion  of  the  spread  jamming  signal  remains 
within  the  frequency  band  entering  the  tracking  loop  with  the  satellite  signal.  Thus  signal  tracking  can  be  stable 
only  if  the  jamming-to-signal  ratio  (JSR)  is  below  the  processing  gain  of  the  spreading  code.  In  general,  a  JSR 
of  greater  than  40dB  is  likely  to  prevent  the  GPS  receiver  from  tracking  the  satellite  signal  and  from  estimating 
its  own  position. 

There  are  situations  where  jammers  may  be  much  stronger  than  the  GPS  signals,  and  are  located  at  close 
physical  proximity  to  the  GPS  receiver.  In  such  a  case,  the  spreading  gain  of  the  spread-spectrum  system  might 
not  be  sufficient  to  decode  the  satellite  data  reliably.  In  fact,  one  of  the  major  difficulties  in  time-space  position 
information  (TSPI)  design  in  defense  applications  is  jamming  rejection.  During  the  project  period,  we  developed 
a  procedure  based  on  the  empirical-mode  decomposition  (EMD)  method  [28],  in  combination  with  the  traditional 
blind  source  separation  (BSS)  technique,  as  a  receiver-based  algorithm  to  suppress  jamming.  For  practical  im¬ 
plementations,  our  algorithm  may  be  integrated  into  the  software  radio  such  as  the  one  described  in  Ref.  [29]. 
Schemes  similar  to  ours  have  been  suggested  for  applications  in  different  contexts,  such  as  biomedical  signal 
processing,  notably  in  Refs.  [30]  and  [31]  where  a  combination  of  EMD  and  independent  component  analysis 
was  employed. 

The  traditional  Fourier  analysis  is  powerful  for  stationary  (or  piecewise  stationary)  signals.  While  the  wavelet 
method  can  handle  nonstationary  signals,  it  is  essentially  an  adjustable  window  Fourier  spectral  analysis  and 
therefore  fundamentally  it  is  also  a  linear  method.  The  method  of  EMD  was  pioneered  by  Huang  et  al.  to 
specifically  deal  with  nonstationary  and/or  nonlinear  signals  [28].  That  jamming  signals  in  GPS  applications 
are  likely  to  be  nonstationary  and  possibly  nonlinear  motivates  us  to  investigate  the  applicability  of  EMD-based 
methods  for  antijamming.  The  block  diagram  of  our  proposed  method  is  shown  in  Fig.  6,  where  the  original 
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Figure  7:  GPS  direct-sequence  spread-spectrum  signal  structure  (C/A  code). 


GPS  data  that  carry  the  satellite  information  are  encoded  and  both  noise  and  jamming  are  present  in  the  received 
GPS  signal.  The  EMD/BSS  combined  procedure  decomposes  the  jammed  GPS  signal  into  a  number  of  modes 
of  distinct  time  scales.  Since  the  code  that  encodes  the  satellite  information  (the  C/A  code)  is  also  known  at  the 
receiver  end,  a  decorrelation  procedure  can  be  used  to  identify  the  modes  that  contain  mostly  the  original  GPS 
data.  We  shall  detail  our  method  and  provide  strong  evidence  that  the  method  works  for  nonstationary  jamming 
of  JSR  of  up  to  45  dB. 

In  the  following  we  provide  a  brief  background  of  the  GPS  signal  structure  and  discuss  several  existing 
antijamming  methods.  We  then  describe  the  EMD  and  BSS  methods  and  their  implementations  and  present 
simulation  results  with  GPS  signals. 

3.2.2  GPS  signal  structure  and  existing  antijamming  methods 

GPS  consists  of  24  satellites  orbiting  at  the  altitude  of  20,183  km,  of  known  positions.  The  signals  transmitted 
by  the  GPS  satellites  are  direct-sequence  spread-spectrum  (DSSS)  signals  that  consist  of  three  portions:  the 
Course  Acquisition  (C/A)  code  on  the  LI  carrier  (1.575GHz),  the  P-code  on  the  LI  carrier,  and  the  P-code  on 
the  L2  carrier  (1.227GHz).  The  C/A  code  has  a  chip-rate  of  1.023MHz  and  a  period  of  1msec,  while  the  P-code, 
for  military  use  only,  has  a  chip-rate  of  10.23MHz  and  a  period  of  1  week.  As  shown  in  Fig  7,  with  Binary 
Phase  Shift  Keying  (BPSK)  modulation,  the  resulting  C/A -code  signal  requires  a  bandwidth  of  2x  1.023MHz  for 
transmission,  and  the  P-code  signal  needs  a  bandwidth  of  2x  10.23MHz. 

A  GPS  receiver  receives  DSSS  signals  from  four  or  more  satellites  and  estimates  the  code-phase  differences 
and/or  carrier-phase  differences  in  order  to  calculate  its  own  position  [2].  As  shown  in  Fig.  8,  the  GPS  receiver 
contains  two  tracking  loops:  code-phase  tracking  loop  and  carrier-phase  tracking  loop.  Code  tracking  and  de¬ 
spreading  are  performed  prior  to  carrier  tracking  since  the  power  of  the  received  GPS  signal  is  much  lower  than 
the  background  noise  and  sufficient  signal-to-noise  ratio  (SNR)  necessary  for  carrier-phase  tracking  can  only  be 
achieved  by  code  despreading.  The  pseudo-random  noise  (PN)  generator  generates  a  PN  sequence  identical  to 
the  C/A  or  the  P  code  to  synchronize  with  the  input  signal  from  a  GPS  satellite.  The  code  correlator  sweeps  the 
uncertainty  ranges  of  the  input-code  phase  at  discrete  steps,  and  detects  the  coarse  synchronization  (acquisition) 
of  the  code  phases  between  the  input  and  local  signals  by  finding  the  maximum  of  the  correlation  function.  The 
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Figure  8:  Block  diagram  of  a  GPS  receiver  that  typically  consists  of  a  code  and  a  carrier-phase  tracking  loop. 


code-phase  tracking  loop  can  track  the  variations  of  the  incoming  code  phase  and  keeps  the  code-phase  alignment 
error  within  an  allowable  limit  after  the  code-phase  acquisition.  From  accurate  tracking  of  the  code  phase,  the 
pseudo-range  time  delay  can  be  obtained  and  the  input  signal  is  despread  to  yield  encoded  navigation  and  timing 
information.  The  despread  signal  is  then  passed  on  to  the  carrier-phase  tracking  loop  for  ranging. 

In  a  hostile  environment  where  jamming  sites  may  be  close  to  GPS  users,  a  huge  JSR  is  possible.  Thus,  in  or¬ 
der  for  GPS  receiver  to  function,  jammer  should  be  either  rejected  before  the  tracking  loops  or  reduced/eliminated 
inside  the  tracking  loops.  For  the  first  case,  the  structure  of  conventional  tracking  loops  does  not  require  change, 
while  for  the  second  case,  novel  algorithms  need  to  be  designed. 

There  are  a  number  of  ways  to  mitigate  the  effects  of  jamming  on  GPS  receivers  before  the  signal  enters  the 
tracking  loops.  These  are  (1)  time-domain  filtering  [32,  33],  (2)  frequency-domain  filtering  [34,  35],  (3)  spatial 
filtering  [24,  36],  and  (4)  time-frequency  filtering  [22,  37].  The  first  two  types  of  filtering  are  conventional.  The 
third  type  typically  uses  adaptive  nulling  antenna,  an  array  of  antenna  elements.  Spatial  filters  have  the  ability 
to  modify  its  reception  pattern,  i.e.,  different  amplification  rates  for  signals  from  different  directions  of  arrival 
(DOA).  Based  on  the  assumption  that  jamming  signals  and  satellite  signals  have  different  DOA,  adaptive  nulling 
antenna  can  emphasize  the  desired  GPS  satellite  signals  and  reject  the  jamming  signals  [24].  The  technique  is 
effective  for  both  narrow-  and  broad-band  jamming  signals.  However,  due  to  multi-path  propagation  of  signals 
and  constraints  on  its  size  and  power,  adaptive  antenna  alone  cannot  provide  an  acceptable  interference  mitiga¬ 
tion.  The  fourth  type  (time-frequency  filtering)  relies  on  the  assumption  that  broad-band  satellite  and  jamming 
signals  have  distinct  time-frequency  signatures.  Once  the  instantaneous  frequencies  of  the  jamming  signals  are 
estimated  from  the  received  signals,  techniques  such  as  time-varying  notch  filter  [22]  or  subspace  projection  [37] 
can  be  used  to  reduce  or  eliminate  the  jamming.  A  possible  drawback  is  that  time-frequency  filters  tend  to  block 
the  frequencies  occupied  by  jamming  signals,  and  thus  also  subtract  the  power  of  satellite  signals  at  those  fre¬ 
quencies  from  the  total  power  of  the  received  signals.  If  jamming  signals  occupy  a  substantial  bandwidth  in  the 
signal  spectrum,  the  filtered  signals  will  be  significantly  distorted.  Thus,  in  general,  this  method  is  more  effective 
for  narrow-band  jamming  signals. 

The  conventional  code-tracking  loop  utilizes  delay-lock  loop  (DLL),  which  offers  little  protection  against 
jamming  and  multipath.  One  way  to  increase  the  robustness  of  DLL  is  to  include  models  (e.g.,  AR  model)  for 
jammer  signals  and  multipath  in  the  DLL.  The  resulting  DLL  can  then  estimate  the  code  delay  using  adaptive 
algorithms,  e.g.,  Kalman  filters  and/or  particle  filters  [38,  39]. 

Our  proposed  EMD/BSS  procedure  is  fundamentally  a  nonlinear  signal-processing  method  and  it  is  therefore 
different  from  these  existing  linear  methods. 
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3.2.3  Nonlinear  antijamming:  proposed  EMD/BSS  methodology 

Our  method  applies  EMD  to  the  received  jammed  GPS  signal  to  decompose  it  into  a  number  of  intrinsic  mode 
functions  (IMFs),  which  can  be  regarded  as  multiple  observations  of  a  random  process.  Assuming  that  these 
observations  are  linear  combinations  of  the  original  signal,  a  BSS  technique  can  be  used  to  extract  the  GPS 
signal  from  the  IMFs. 

EMD  method.  Traditional  methods  such  as  the  Fourier  spectral  analysis  assume  stationarity  anVi/or  approxi¬ 
mate  the  physical  phenomena  with  linear  models.  These  approximations  may  lead  to  spurious  components  in  the 
time-frequency  distribution  diagrams  if  the  underlying  signal  is  nonstationary  and  nonlinear.  Empirical  Mode  De¬ 
composition  (EMD)  is  a  technique  [28]  to  deal  specifically  with  nonstationary  and  nonlinear  signals.  Given  such 
a  signal,  the  method  adaptively  decomposes  it  into  a  number  of  modes  (IMFs)  that  are  topologically  equivalent  to 
amplitude-  and  frequency-modulated,  sinusoidal  signals.  In  the  analytic-signal  representation,  the  modes  corre¬ 
spond  to  proper  rotations  [28].  Thus  the  EMD  method  naturally  yields  estimates  of  the  significant  instantaneous 
frequencies  embedded  in  the  signal,  by  performing  the  Hilbert  transform  on  each  IMF.  The  ease  and  accuracy 
with  which  the  EMD  method  processes  nonstationary  and  nonlinear  signals  have  led  to  its  widespread  use  in 
various  applications  such  as  seismic  data  analysis  [28],  chaotic  time  series  analysis  [40,  41],  and  meteorological 
data  analysis  [42],  etc. 

Given  a  signal  X(t)y  the  EMD  method  focuses  on  the  level  of  local  oscillations  and  decomposes  the  signal 
into  a  finite  and  often  a  small  number  of  fundamental  oscillatory  modes.  The  bases  (IMFs)  into  which  the  signal 
is  decomposed  are  obtained  from  the  signal  itself,  and  they  are  defined  in  the  time  domain.  They  are  of  the 
same  length  as  the  original  signal  and  preserve  the  frequency  variations  with  time.  The  base  modes  can  be 
made  approximately  complete  and  nearly  orthogonal  with  respect  to  each  other.  Here,  completeness  implies 
that  the  original  signal  can  be  reconstructed  without  any  loss  of  information  by  simply  summing  up  the  IMFs. 
Thus,  the  IMFs  can  be  viewed  as  linear  components  of  the  original  or  source  signal  X(t).  In  order  to  achieve 
this,  two  conditions  need  to  be  satisfied:  (1)  the  total  number  of  extrema  of  IMF(£)  be  equal  to  the  number 
of  zero  crossings,  and  (2)  the  mean  of  the  upper  envelope  IMFU(£)  and  the  lower  envelope  IMF  ft)  be  zero. 
The  process  to  obtain  the  IMFs  from  the  signal  is  called  sifting ,  which  consists  of  the  following  steps:  (1) 
identification  of  the  extrema  of  X(t ),  (2)  interpolation  of  the  set  of  maximal  and  minimal  points  (say,  by  using 
cubic  splines)  to  yield  an  upper  envelope  Xu(t)  and  a  lower  envelope  Xft ),  respectively,  and  their  average 
m(t)  =  [Xu(t)  4-  Xft)}/ 2,  (3)  subtraction  of  the  average  from  the  original  to  yield  d(t)  =  X(t)  -  m(£),  and 
(4)  repetition  of  steps  (1-3)  until  d(t)  satisfies  the  two  conditions  for  being  an  IMF.  Once  an  IMF  is  generated, 
the  residual  signal  r(t)  =  X(t)  -  IMF](£)  is  regarded  as  the  original  signal,  and  steps  (1-4)  are  repeated  to  yield 
the  second  IMF,  and  so  on.  The  procedure  is  complete  when  either  the  residual  function  becomes  monotonic, 
or  when  the  amplitude  of  the  residue  falls  below  a  pre-determined  small  value  so  that  further  sifting  would  not 
yield  any  useful  components.  These  features  guarantee  the  computation  of  a  finite  number  of  IMFs  within  a  finite 
number  of  iterations.  The  outcome  of  the  EMD  procedure  is  the  following  decomposition  of  the  original  signal: 

n 

X(()  =  JjMFt(<)  +  r(0  '  (24) 

1=1 

where  IMF*(£)  is  the  ith  IMF,  n  is  the  total  number  of  IMFs,  and  r(t)  is  the  final  residue  that  has  near  zero 
amplitude  and  frequency. 

Figure  9  shows  the  IMFs  obtained  from  a  simulated  GPS  signal  to  which  a  zero-mean  and  unit-variance 
Gaussian  noise  and  a  stationary,  sinusoidal  jamming  signal  are  added.  The  GPS  signal  has  the  amplitude  of  2  and 
the  jamming  amplitude  is  200  (corresponding  to  JSR  of  40dB).  The  top  left  panel  shows  the  noisy  and  jammed 
GPS  signal,  and  the  seven  remaining  panels  show  the  seven  significant  IMFs.  We  note  that  the  first  IMF  in  the  top 
right  panel  visually  resembles  the  original  GPS  signal.  This  example  thus  illustrates  that  for  stationary  jamming, 
the  EMD  method  is  capable  of  directly  separating  the  GPS  signal  from  the  jamming.  However,  this  appears  not 
to  be  the  case  for  nonstationary  jamming,  as  shown  in  Fig.  10,  where  the  top  left  panel  is  the  noisy  and  jammed 
GPS  signal.  The  jamming  is  a  frequency-modulated  signal  with  normalized  frequency  increased  linearly  from  0 
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Figure  9:  IMFs  of  a  noisy  and  jammed  GPS  signal,  where  the  jamming  is  stationary  and  the  noise  is  Gaussian. 
The  first  IMF  (upper  right  panel)  is  approximately  the  original  GPS  signal. 


to  0.5  in  the  time  interval  considered.  We  find  that  none  of  the  seven  IMFs  shown  (in  the  remaining  7  panels) 
looks  like  the  original  GPS  signal,  which  is  in  fact  embedded  in  all  IMFs.  It  is  thus  necessary  to  invoke  a  proper 
BSS  procedure  to  extract  the  GPS  signal. 

Blind  Source  Separation.  Blind  Source  Separation  (see  [43]  for  a  review)  is  a  method  to  extract  basic  source 
signals  from  several  observed  mixtures.  It  is  considered  “blind”  because  the  source  signals  are  not  observed 
and  no  information  is  available  about  the  observed  mixtures.  An  a  priori  assumption  is  that  the  sources  are 
independent  of  each  other  and  they  are  not  all  white  noise.  BSS  is  particularly  suitable  for  situations  where 
different  observations  of  the  same  sources  are  received  from  different  sensors.  The  BSS  technique  has  found 
applications  in  areas  such  as  communications  and  biomedical  signal  processing  [44,  45,  46]. 

Consider  the  time  varying  vector  of  observations 

X(t)  =  [Xl(t),..,Xk(t)]T  (25) 


obtained  as  a  mixture  from  k  sources 

S(t)  =  [Si  (t), ...  Sk(t)]T  and X(t)  =  AS(i)  (26) 

where  A  is  the  mixing  matrix.  We  reconstruct  the  sources 

Y  (t)  =  B  X(t)  (27) 

where  B  is  a  matrix  to  be  determined,  by  adopting  the  criterion  of  minimizing  mutual  information.  This  requires 
the  estimation  of  joint  entropy,  which  is  computationally  expensive,  and  therefore  we  adopt  the  idea  of  Gaussian 
Mutual  information  [47]  that  requires  computation  only  up  to  second-order  characteristics. 
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Figure  10:  IMFs  of  a  noisy  and  jammed  GPS  signal,  where  the  jamming  is  nonstationary  and  the  noise  is  Gaus¬ 
sian.  In  this  case,  the  original  GPS  signal  is  spread  over  all  IMFs. 


3.2.4  Simulation  Results 

For  simulation  purposes,  we  have  used  the  following  baseband  model  for  received  GPS  signals: 

r(t)  =  c(t)d(t)  -F  j(t)  +  n(t),  (28) 

where  r(t)  is  the  received  signal,  c(t)  is  the  spreading  sequence,  d(t)  is  the  transmitted  GPS  information  symbol, 
n(t)  is  the  noise,  and  j(t)  is  the  jamming.  In  the  simulation,  the  length  of  spread  sequence  is  varied  from  10  to 
20  bits.  For  stationary  jamming  simulations,  pure  sinusoids  of  frequency  0.2  to  0.5  times  the  sampling  frequency 
are  chosen.  Nonstationary  jamming  is  modeled  as  a  chirp  signal  with  frequency  changing  from  0  to  0.5  times 
the  sampling  frequency.  The  amplitude  of  the  transmitted  GPS  signal  is  normalized  to  unity,  while  the  jamming 
amplitude  is  varied  so  that  the  maximal  JSR  is  45dB. 

The  code  for  EMD  is  adapted  from  the  one  developed  by  Rilling  et  al.  [Matlab  codes  (G.  Rilling,  P.  Flandrin 
and  P.  Goncalves),  http://perso.ens-lyon.fr/patrick.flandrin/emd.html].  During  a  run,  typically  between  8  and  10 
IMFs  are  generated.  As  can  be  seen  from  the  time-series  plots  in  Fig.  9,  from  the  IMFs  of  stationary  jammed  GPS 
signal  we  can  visually  distinguish  the  jammer  from  the  GPS  signal.  The  GPS  signal  is  captured  almost  entirely 
in  the  first  mode,  whereas  the  jammer  is  captured  in  the  2nd  and  3rd  modes.  In  this  case,  the  BSS  procedure  is 
not  necessary.  To  recover  the  GPS  signal,  it  is  a  relatively  simple  matter  of  correlating  the  mode  that  contains 
the  GPS  signal  with  the  PRN  code.  For  nonstationary  jammer,  as  in  Fig.  10,  both  the  GPS  and  the  jamming 
are  mixed  in  all  the  IMFs.  In  this  case,  it  is  necessary  to  use  the  BSS  to  extract  the  GPS  signal.  In  particular, 
the  IMFs  can  be  regarded  as  multiple  observations  of  the  received  signal.  We  adapted  the  code  from  the  BLISS 
project  [BSS  Demonstration  Code  (http://www-lmc.imag.fr/SMS/SASI/bliss.html)]  for  BSS. 

Figure  1 1  shows  the  result  of  extracting  the  GPS  signal  in  the  presence  of  stationary  jamming,  where  the  two 
top  panels  show  the  original  GPS  and  the  original  jamming  signals,  respectively,  and  the  two  lower  panels  show 
the  extracted  GPS  and  jamming  signals.  A  typical  result  with  nonstationary  jamming  is  shown  in  Fig.  12.  It  is 
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Figure  11:  Result  with  stationary  jammer.  Upper  panels:  original  GPS  signal  and  jammer.  Lower  panels: 
extracted  GPS  signal  and  jammer. 


apparent  that  our  combined  EMD/BSS  methodology  is  capable  of  estimating  the  GPS  signal  in  the  presence  of 
strong  jamming,  stationary  or  nonstationary. 

To  quantify  the  “goodness”  of  our  method,  we  examine  the  bit  error  rate.  In  particular,  information  bits 
associated  with  the  extracted  GPS  signal  are  compared  with  those  in  the  original  GPS  signal.  A  typical  result  for 
one  C/A  code  is  shown  in  Fig.  13,  where  the  ratio  of  bits  received  correctly  to  the  number  of  bits  transmitted 
versus  the  jamming  signal  amplitude  is  plotted.  The  upper  and  lower  traces  show  the  ratio  for  stationary  and 
nonstationary  jamming,  respectively.  We  see  that  for  stationary  jamming,  the  information  bit  can  be  extracted 
with  certainty.  For  nonstationary  jamming  of  JSR  of  as  high  as  45dB,  the  percentage  of  correctly  extracted  bit  is 
above  80%.  Since  one  GPS  bit  (0  or  1)  is  modulated  using  20  C/A  codes,  setting  a  conservative  threshold,  say 
at  unity,  for  the  ratio  between  the  number  of  correctly  estimated  bits  and  that  of  the  wrong  bit  can  guarantee  the 
extraction  of  the  correct  bit.  For  example,  suppose  the  original  bit  is  1  and  the  JSR  is  40dB,  then  out  of  the  20 
C/A  codes,  1  will  appear  approximately  16  times  and  0  (the  wrong  bit)  will  appear  about  4  times.  The  ratio  is 
then  4,  which  can  be  distinguished  from  1  almost  certainly.  Figure  13  suggests  that  our  EMD/BSS  method  can 
perform  to  extract  GPS  satellite  information  in  the  presence  of  nonstationary  jamming  of  JSR  up  to  about  45dB. 

3.2.5  Discussions 

Intentional  jamming  of  GPS  signals  is  a  serious  concern  especially  for  military  applications.  The  nature  of 
the  GPS  signal,  i.e.,  its  extreme  low  signal  intensity,  makes  it  vulnerable  to  jamming.  In  applications  where 
dependence  on  GPS  is  high,  a  jammed  GPS  signaPcould  have  disastrous  consequences.  Thus  it  is  of  tremendous 
interest  and  importance  to  explore  practical  antijamming  schemes. 

We  proposed  and  tested  a  simple  yet  effective  antijamming  algorithm  which  can  be  implemented  in  the 
receiver  stage  after  the  acquisition  of  the  GPS  signal.  The  algorithm  makes  use  of  two  available  techniques  for 
data  analysis:  the  Empirical  Mode  Decomposition  method  and  the  Blind  Source  Separation  method  based  on 
Gaussian  Mutual  Information.  We  showed  that  the  algorithm  is  capable  of  accurate  estimation  of  GPS  signal  bits 
in  the  presence  of  stationary  or  nonstationary  jammers.  In  particular,  simulation  results  indicate  that  the  GPS 
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Figure  12:  Result  with  nonstationary  jammer.  Upper  panels:  original  GPS  signal  and  jammer.  Lower  panels: 
extracted  GPS  signal  and  jammer. 
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Figure  13:  For  one  C/A  code,  the  percentage  of  correctly  extracted  number  of  GPS  signal  bits  versus  the  jamming 
amplitude.  The  upper  trace  is  for  stationary  jammer  and  the  lower  for  nonstationary  jammer.  Since  one  GPS 
information  bit  is  encoded  using  20  C/A  codes,  for  JSR  of  up  to  45dB,  the  correct  bit  can  be  obtained  with 
certainty. 
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data  can  be  decoded  accurately  with  no  error  for  nonstationary  jammer  of  JSR  of  up  to  about  45dB. 

For  nonstationary  jammer  of  JSR  above  50dB,  the  performance  of  our  method  deteriorates.  The  solution 
may  be  to  combine  this  algorithm  with  spatial  filtering  methods  such  as  those  described  in  [24,  36].  Another 
possibility  is  to  integrate  Kalman  or  particle  filtering  in  the  algorithm  for  estimating  the  phase/code  delay.  The 
model  for  such  a  system  would  assume  the  presence  of  jamming  interference  and,  when  combined  with  our 
algorithm,  might  enable  GPS  signal  processing  in  the  presence  of  exceedingly  strong  jammer. 

3.3  Precise  GPS  positioning  of  moving  objects  in  noisy  environment  by  particle  filters 

3.3.1  Background 

The  problem  of  estimating  the  time-varying  state  of  a  system  based  on  experimental  measurements  or  obser¬ 
vations  finds  many  applications  in  physics,  biology,  and  engineering.  Examples  include  quantum  state  recon¬ 
struction  and  purity  estimation  [48,  49],  noise  reduction  and  state  reconstruction  in  chaotic  dynamical  systems 
[50,  51],  estimation  of  bounds  on  isocurvature  perturbations  in  the  Universe  and  on  cosmic  strings  from  cos¬ 
mic  microwave  background  and  large-scale  structure  data  [52,  53],  gravitational  wave  signal  analysis  [54,  55], 
macromolecular  structure  determination  [56,  57],  prediction  of  protein-protein  interactions  from  genomic  data 
[58],  tracking  and  positioning  problems  [59],  etc.  In  general,  a  system  model  that  evolves  the  state  vector  (x) 
in  time  is  needed,  so  is  an  observational  model  that  relates  an  observation  vector  ( y )  to  the  state  vector.  In  any 
realistic  application  both  noise  and  model  uncertainties  exist,  rendering  necessary  a  probabilistic  treatment  of  the 
estimation  problem.  That  is,  one  can  evolve  x  according  to  the  system  model,  and  make  corrections  to  or  update  x 
based  on  the  available  y.  The  quantity  of  interest  is  the  posterior  probability  density  function  (pdf)  p{x\y),  given 
all  available  observations  y.  The  standard  approach  to  addressing  this  problem  is  Bayesian  inference  [60],  which 
leads  to  the  classical  Kalman  filter  when  both  the  system  and  the  observational  models  are  linear.  For  nonlinear 
problems,  a  viable  approach  is  sequential  Monte-Carlo  simulation  (or  particle  filter)  [61,  62,  63],  which  uses  a 
set  of  random  samples  to  approximate  the  posterior  pdf  p(x\y).  The  approximated  pdf  evolves  and  is  corrected 
by  the  observation  based  on  the  Bayes  rule.  If  the  number  of  samples  is  sufficiently  large,  the  approximation 
approaches  the  optimal  Bayesian  estimate.  Due  to  the  constant  improvement  of  modem  computing  technol¬ 
ogy,  the  sequential  Monte-Carlo  approach  has  begun  to  find  significant  applications  in  science  and  engineering 
[62,  63,  64]. 

A  fundamental  question  in  sequential  Monte-Carlo  simulations  is  how  the  precision  of  the  estimated  state 
vector  depends  on  noise  in  the  system  and  in  the  observation.  Another  issue  of  significant  practical  interest  is 
how  to  deal  with  occasional  but  large  disturbances,  or  outliers ,  in  the  observation.  During  the  project  period,  we 
addressed  these  two  related  problems.  In  particular,  we  derived  and  verified  a  self-consistent  equation  that  relates 
the  covariance  matrix  of  the  samples,  which  determines  the  precision  of  the  state  estimate,  to  the  covariance 
matrices  of  the  system  noise  and  of  the  observational  noise.  In  addition,  we  proposed  a  robust  sequential  Monte- 
Carlo  scheme  to  overcome  the  effect  of  outliers.  In  this  regard,  previous  approach  includes  using  heavy-tailed 
error  distribution  to  improve  the  state-space  models  so  that  they  react  quite  flexibly  to  changes  in  points  or  edges, 
but  still  provide  smooth  fits  in  other  regions  [65].  Leave-k-out  diagnostic  is  used  to  detect  a  series  of  consecutive 
outliers  for  a  linear  state  space  model  [66].  It  uses  all  residual  observations  in  the  time  span  to  check  whether 
a  series  of  consecutive  observations  are  jointly  outlying,  thus  it  is  actually  “off-line”.  Our  idea  is  to  detect  the 
outliers  from  the  previous  knowledge  about  the  system  and  then  to  eliminate  them  in  the  sequential  Monte-Carlo 
implementation.  Simulations  using  a  precise  GPS  positioning  problem  demonstrate  the  power  of  our  scheme. 
We  expect  our  results  to  have  significant  impacts  on  problems  where  the  underlying  system  and/or  experimental 
observations  are  subject  to  outliers.  For  example,  the  observation  of  a  star  or  a  galaxy  may  be  corrupted  by  the 
drastic  activity  of  another  celestial  body  in  a  short  period.  In  biological  physics,  macromolecular  structure  is 
inferred  indirectly  from  various  measurements,  e.g.  nuclear  magnetic  resonance  spectra,  X-ray  reflections,  or 
homology-derived  restraints,  which  can  easily  contain  outliers  [56].  In  GPS  (global  positioning  system)-based 
precise  positioning  problems,  GPS  signals  may  be  disturbed  by  sudden  and  large  jamming. 
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In  the  following,  we  outline  the  basic  steps  of  sequential  Monte-Carlo  method  and  derive  a  self-consistent 
equation  governing  the  dependence  of  estimation  error  on  noise.  We  then  present  a  robust  sequential  Monte-Carlo 
scheme  for  mitigating  the  effect  of  measurement  or  observational  outliers  and  provide  numerical  support. 

3.3.2  Sequential  Monte-Carlo  method 

Let  y( 0  :  t)  =  {*/(£'),£'  =  *o(=  0),  t\,  t<i,  •  •  •  ,  tk{=  0}  be  the  observations  from  time  0  to  time  t ,  which  are  not 
necessarily  equidistant  in  time.  We  seek  to  obtgin  the  posterior  pdf  p[x(t)\y(0  :  t)].  The  state  and  observational 
equations  are 


x(t)  =  /[x(ifc_i),w(f/k_i)l,  (29) 

y(t)  =  ff[x(t)»e(0]»  (30) 


which  describes  the  evolution  of  the  state  and  maps  the  state  to  the  observational  vector,  respectively,  /  and  g 
can  be  nonlinear  functions.  The  processes  v(t)  and  e(t)  represent  random  fluctuations  (e.g.,  noise,  uncertainties, 
outliers,  etc.)  in  the  system  and  in  the  observation,  respectively.  Often,  in  an  application  the  distribution  of  the 
initial  state  can  be  obtained  by  considering  the  specific  physics  involved.  It  is  thus  reasonable  to  assume  that 
this  distribution  is  available.  The  pdf  p[x(t)\y(0  :  £)]  can  then  be  obtained  recursively  by  prediction  through  the 
dynamical  equation  (29)  and  likelihood  correction  through  the  observational  equation  (30).  In  particular,  given 
the  pdf  p[x(tk-\)\y(0  :  4-i)]  at  time  t*_j,  the  prediction  step  uses  the  dynamic  equation  (29)  to  obtain  the  prior 
pdf  of  the  state  at  time  t  via  the  Chapman-Kolmogorov  equation 

p[x(t)\y(0  :  =  J  dx(tk- 1)  •  p[z(t)|x(t*_i)j  •  p[x(*jt_i)|y(0  :  k-i)],  (31) 


where  p[x(t)\x(tk-i),y(0  :  tk-i)]  =  p[x(£)|x(ffc_i)]  is  used.  At  time  t ,  a  new  measurement  y(t)  becomes 
available,  which  can  be  used  to  correct  the  prior  pdf  via  the  Bayes  rule 


p[x(<)|y(0  :  t)}  = 


p[y(0b(0]pb(0ly(Q :  4-i)] 

p[y(Oly(o :  **-i)] 


(32) 


where 

p[y(t)  b(o :  4-i)l  =  J  p[y(0l*(01  pMOIy(o :  tk_x))dx 

depends  on  the  likelihood  function  p[y(t)\x(t)].  The  recurrence  relations  (31)  and  (32)  form  the  basis  for  optimal 
Bayesian  solution.  For  Gaussian  noise,  when  /  and  g  are  linear  functions,  the  recurrence  relation  can  be  solved 
analytically,  which  is  the  classical  Kalman  filter.  For  nonlinear  functions  /  and  g ,  a  linearization  technique  is 
viable  which  leads  to  the  so  called  extended  Kalman  filter  [67].  Unscented  Kalman  filter  deliberately  selects  a 
set  of  points  and  propagate  them  through  the  nonlinearity  to  estimate  the  Gaussian  approximation  [68].  While 
for  more  general  cases,  the  approach  of  sequential  Monte-Carlo  simulation  is  desirable  [61, 62,  63,  64]. 

Let  {xj(f),  denote  a  random  measure  that  characterizes  the  posterior  pdf  p[x(t)\y(0  :  £)],  where 

{xj(£),  i=l,  ...  ,N}  is  a  set  of  support  points  with  associated  weights  {iut(f),  i  =  1,  ...  ,7V}.  The  weights  are 
normalized  that  Wi(t)  =  1.  The  posterior  pdf  can  be  approximated  as 


N 

p[x(t)|y(0  :  t)]  ~  ^ ~2wi(t)S{x{t )  -  Xi(t)), 
1=1 


where  6(x)  is  the  Dirac  delta  function  such  that  S(x)  =  0  if  x  ^  0  and  fi(x)dx  =  1  if  (xj,  X2)  contains  0. 
The  average  of  an  arbitrary  function  f(x)  can  be  simplified  as 

(/[*(0)>  =  j  f(x)p[x\y( 0  :  t)]dx  =  ^  w«(0/[*i(01- 
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The  weights  are  chosen  using  the  principle  of  importance  sampling  [69].  In  particular,  given  an  arbitrary  pdf 
p{x),  it  may  be  difficult  to  draw  samples.  Suppose  for  an  alternative  pdf  q(x ),  samples  can  be  drawn  relatively 
easily.  Letting  X{  ~  q(x)  (z  =  1, ...,  N)  be  samples  drawn  from  some  <?(•),  the  importance  density ,  we  obtain  the 
following  weighted  approximation: 

N 

p(x)  ~  ^  W{S(x  —  Xi), 

1=1 

where  Wj  oc  p(xi)/q(xi)  is  the  normalized  weight  of  the  z’th  sample. 

Now  consider  the  joint  probability  p[x( 0  :  |?/(0  :  £)].  In  case  of  independent  noise  samples,  we  can  write 

k 

p[z(0  :  f)|y(0  :  f)]  oc  p[x(0)|y(0)]  p[y((_;)|x(^)]p[x(<i)|x((j_1)]. 


Thus 


p[x(0  :  <)|y(0  :  t)}  =  p[x(t),x(0  :  tk_i)\y(t),y(0  :  tfc-i)] 

p[y(«)|x(t)]p[x(t)|x(4-i)]. 


p[y(<)|y(0  :  <fc-i)] 


■p[x(0  :  i ) |y(0  :  <fc-i)], 


Assume  the  posterior  distribution  p[x(0  :  ) |y(0  :  tfc-i)]  is  approximated  by  {x<(0  :  tk-\),  Wi{tk_ i)}£Lj, 

given  a  new  observation  y(f),  the  objective  is  to  obtain  an  approximation  {x<(0  :  t),  iu»(£)}£Li  for  p[x(0  : 
f)|y(0  :  <)],  such  that  the  estimation  of  quantities  of  interest  at  time  t  can  be  calculated.  The  sequential  Monte 
Carlo  scheme  is  to  generate  a  sample  Xi(t)  and  append  it  to  Xi(0  :  t^-i)  to  form  Xj(0  :  t ),  and  update  the  weight 
Wi(tk- 1)  to  Wi{t). 

If  the  importance  function  y[x(0  :  i )  |y  (0  :  t )]  can  be  factorized  as 


y[x(0  :  t)|y(0  :  t)}  =  y[x(f)|x(0  :  ffc-i),y(0  :  £)]  *  <7[x(0  :  ffc_i)|y(0  :  (33) 


and  x,(<)  is  sampled  from  y[x(f)|xj(0  :  <*-1)1  y(0  :  £)],  the  weight  of  the  trajectory  x<(0  :  t)  is 

(t)  P[x»(0  :  <)|y(0  :  t) |  =  P[y(0|xj(t)]p[xi(f)|xf(4_i)) _  x  p[x,(0  :  tfc-i)ly(0  :  £fc_  1)] 

^  ?[xi(0  :  f)|y(0  :  £)]  y[xi(<)|xi(0  :«*_]), y(0  :  f)]p[y(f)|y(0  :  4-i)]  ?[xj(0  :  <fc_,)|y(0  : 

p[y(t)|xj(<)]p[xi(^)|xi(ffc_1)] 

^  y[xi(<)|xi(0  :  y(0  :  £)]  Wl  k~1 

where  p[y{t)\y(0  :  tk-i)]  is  omitted  since  it  is  common  to  all  samples.  A  convenient  choice  for  the  importance 
density  is  the  following  prior  importance  density  [63,64]:  ^[x(<)|xi(0  :  ^_i),y(0  :  £)]  =  p[x(t)\xi(tk-\ )],  with 
which  the  weight  updating  equation  becomes 


Wi{t)  oc  wi(<fc_i)p[y(f)|xi(f)],  (34) 

and  the  posterior  filtered  density  p[x(t)\y(0  :  £)]  can  be  approximated  as 

N 

p[x(£)|y(0  :  <)]  ~  w«(t)<)[x(<)  -  Xj(f )].  (35) 

»=l 

From  a  numerical  point  of  view,  the  above  analysis  can  be  implemented  as  follows.  First  generate  N  samples 
{xj(0)|z  =  1, •  •  •  ,  N}  from  the  distribution  of  y(0)  as  given  in  Eq.  (30).  Each  sample  has  a  weight  of  1/7V. 
Each  sample  x;(0)  then  evolves  according  to  the  dynamical  equation  (29)  by  considering  the  noise  v  to  get  the 
value  at  time  t\,  e.g.  Xi(t\ ),  and  the  weights  are  updated  via  Eq.  (34).  The  estimation  of  the  state  at  time 
t\  is  (x(t]))  =  YliL\  Wi(£i)xi(£i).  During  the  evolution,  it  may  occur  that  there  are  disproportionally  fewer 
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samples  about  x*  than  determined  by  weight  wx.  To  avoid  this,  a  sample  importance  resampling  procedure  [63] 
can  be  applied.  That  is,  we  can  generate  a  new  set  of  samples  {x*,  w *  }^Ll  from  the  samples  {x*,  with 

probability  being  their  weights,  i.e.  each  time,  the  probability  to  draw  sample  xt  is  its  weight  wx  (note  that  the 
weights  are  normalized  that  wi  =  !)•  The  weights  w *  for  the  new  samples  are  then  set  as  1/N.  As  a  result  of 
this  resampling  procedure,  the  weight  Wi  of  a  sample  is  represented  by  the  number  of  duplications  of  the  sample, 
thus  the  statistics  of  the  samples,  e.g.  mean  value,  covariance,  etc.,  are  unchanged  in  the  large  N  limit.  The 
resampling  step  automatically  concentrates  the  samples  in  regions  of  interest  and  effectively  discards  samples 
with  low  weight.  However,  this  may  result  in  overlaps  for  some  samples.  For  example,  if  one  sample  has  a 
very  large  weight,  after  resampling,  it  may  have  many  duplications,  which  leads  to  a  degeneration  problem.  To 
overcome  this  difficulty,  a  regularization  process  can  be  applied:  A  small  random  vector  is  added  to  each  sample 
as  a  perturbation:  xx  <—  x*  -f  hDtx,  where  ex  follows  the  standard  normal  distribution,  D  is  such  that  DDT  =  5 
( D 1  is  the  transpose  of  D).  The  matrix  S  is  the  empirical  covariance  matrix  of  the  samples  before  resampling, 
and  the  quantity  h  is  a  regularization  parameter  [63,  70].  The  samples  again  propagate  via  the  dynamical  equation 
(29)  to  yield  the  values  at  next  time  step.  The  process  continues  until  a  desired  time  span  for  estimation  is  reached. 

3.3.3  Noise  dependence  of  estimation  error 

When  noise  of  the  system  is  stationary,  i.e.  the  covariance  matrices  Ev  and  Ee  for  the  process  noises  v  and  e 
in  Eqs.  (29)  and  (30)  are  constant  in  time,  the  samples  can  evolve  into  a  “steady”  state  and  their  covariance 
matrix  can  be  obtained,  which  is  proportional  to  the  estimation  error.  Suppose  at  time  tk- 1  the  covariance  matrix 
of  the  samples  is  Ex(4_i),  which  is  unknown.  Since  the  dynamical  equation  /  is  known,  after  propagating 
through  Eq.  (29),  the  covariance  matrix  Ey  of  the  samples  at  time  4  can  be  expressed  in  terms  of  Ex(4-i)  and 
Ev,  which  reads  E/(Ex(^_i),  Ev).  To  make  use  of  the  correction  step  [Eq.  (32)],  we  solve  x  from  Eq.  (30): 
x  =  g~l{y,e).  Therefore,  for  a  given  observation  y(t)>  the  covariance  matrix  for  x(£),  from  the  observational 
point  of  view,  can  be  obtained  as  E s[y(t),  Ee].  Usually,  E5  depends  mainly  on  Ee,  and  has  little  dependence  on 
y(t ),  thus  E5  is  merely  constant  in  time  and  can  be  calculated  using  the  initial  observation  y(0).  The  correction 
procedure  is  equivalent  to  a  modulation  posted  by  a  distribution  with  covariance  matrix  E5  on  a  distribution  with 
covariance  matrix  E y.  Suppose  both  distributions  are  Gaussian,  the  resulting  distribution  is  also  Gaussian  but 
with  covariance  matrix  {E J1  -F  E^1}-1.  The  resampling  step  does  not  change  the  covariance  matrix,  and  the 
regularization  step  simply  adds  a  factor  of  1  -f  h 2 .  Thus  we  have  E x(t)  =  {[EJ1  4-  Ej1}-1^  4-  h2).  In  the 
steady  state,  we  have  E x{t)  ~  Ex(£*.-i),  leading  to  the  following  self-consistent  equation: 

EJ1  +  E71  =  (1  +  h2)Z,-\  (36) 

which  determines  the  covariance  matrix  of  the  samples,  or  the  posterior  pdf  p[x(t)\y(0  :  £)],  for  given  dynamical 
and  observational  noise  levels.  Note  that  Ey  is  a  function  of  Ex  and  Ev.  For  certain  cases,  E y  can  be  expressed 
explicitly  in  terms  of  Ex  and  Ev,  which  can  be  used  to  further  simplify  the  above  equation.  For  example,  for 
linear  dynamical  systems,  /  =  y/ax  4-  Vbv,  Ey  =  aEx  4-  6EV,  Eq.  (36)  becomes: 

(aE*  +  6E„)-1  +  Ej1  =  (1  +  /i2)E“\  (37) 

From  Eq.  (36),  the  dependence  of  Ex  on  Ev  and  Ee  can  be  obtained,  which  can  be  used  to  identify  the  “leading 
term”  of  the  noise  source,  i.e.,  which  noise  term  has  the  most  influence  on  Ex,  and  therefore  on  the  estimation 
precision.  This  information  can  be  useful  for  improving  the  estimation  precision  by  suppressing  the  leading  noise 
source.  In  practice,  due  to  the  nonlinearity  of  function  gy  an  explicit  expression  of  E5  is  not  always  possible,  and 
a  Monte  Carlo  scheme  is  viable,  i.e.,  draw  a  set  of  samples  {ej}  from  the  distribution  of  e  [Eq.  (30)];  for  each 
sample  e*,  calculate  xx  as  xt*  =  <7-1(y(0),ei),  then  E s  can  be  approximated  by  the  covariance  matrix  of  the 
samples  {xi}. 
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3.3.4  Robust  sequential  Monte-Carlo  scheme  for  mitigating  outliers 

The  above  scheme  of  sequential  Monte-Carlo  simulation  works  well  for  stationary  noise.  In  the  presence  of 
nonstationary  disturbances,  e.g.,  outliers,  the  weight-updating  scheme  needs  to  be  improved.  To  be  concrete, 
we  treat  outliers  in  the  observation,  which  can  lead  to  a  larger  covariance  matrix  £5.  This  will  cause  a  larger 
estimation  error  [Eq.  (36)].  Thus,  if  the  outliers  can  be  detected  and  then  discarded,  the  elements  of  the  covariance 
matrix  of  the  observations  can  be  reduced.  Our  idea  is  to  first  calculate  the  empirical  distribution  of  the  Monte- 
Carlo  samples,  and  then  compare  the  observation  with  this  distribution.  If  the  prediction  of  the  observation  is 
close  to  the  center  of  the  samples,  the  observation  is  likely  to  be  true  and  it  should  be  accounted  for  in  the 
estimation  of  the  state.  However,  if  the  prediction  deviates  from  the  center  of  the  samples,  it  is  less  reliable 
and  should  therefore  be  counted  less  [71].  Quantitatively,  it  is  convenient  to  introduce  a  contribution  factor  a  to 
characterize  this  effect.  Let  Wi{t)  =  Wi(tk-i)p[y(t)\xi(t)\.  We  modify  Eq.  (34)  as 

wi(t)  =  (1  -  +  avH(t)/^2Sj(t).  (38) 

j 

Generally,  the  optimal  value  of  a  depends  on  the  distribution  of  the  samples  and  on  the  prediction  of  the  ob¬ 
servation  in  a  sophisticated  way.  Our  strategy  for  choosing  a  is  the  following.  After  propagating  the  samples 
through  the  dynamical  equation,  we  calculate  Ex.  The  covariance  matrix  £v  of  the  dynamical  noise  and  £e  of  the 
observational  noise  are  assumed  to  be  known.  Define  ft  =  y/Tr(£v)/Tr(£e),  Ax  =  ( xis  -  (x))/\/7>(£x), 
where  (x)  is  the  average  of  the  samples,  and  xis  is  the  least-squares  estimation  of  the  state,  which  minimizes 
the  square  error  of  the  observations  (this  is  the  case  when  the  number  of  observations  is  more  than  the  number  of 
unknowns)  [72].  We  propose  the  following  criterion  for  choosing  a: 

1,  l|Ax||  <  co, 

I  U-«ife(££S,)J  +  A  <»  <  »Ai||  <  c, 

*  <  »A*II  s  <*. 

,0,  ci  <  ||Ax||, 

where  Co  ~  1,  cj  ~  3.5,  C2  ~  7,  and  the  optimal  values  can  vary  for  different  situations.  Note  that  ||Ax||  is 
the  distance  between  the  estimation  x^s  obtained  from  the  observations  and  the  mean  value  of  the  samples  (x), 
normalized  by  the  “standard  deviation.”  If  ||  Ax||  <  1,  the  estimation  is  within  the  range  of  the  standard  deviation 
and  is  reliable.  If  ||  Ax||  is  in  the  range  of  one  standard  deviation  to  three  standard  deviations,  the  observation  is 
less  reliable.  Since  there  is  noise  in  the  dynamics,  the  samples  may  themselves  contain  some  error.  The  quantity 
(i  is  introduced  to  account  for  such  uncertainties.  If  ||  Ax||  is  even  larger,  the  weight  for  the  observation  decreases, 
and  at  a  certain  point,  say,  beyond  seven  standard  deviations,  the  observation  is  deemed  as  outliers  and  the  weight 
a  is  set  to  be  zero. 

3.3.5  Numerical  support 

To  substantiate  our  ideas,  we  considered  a  synthetic  two-dimensional  GPS  positioning  problem  of  moving  object, 
say  a  car,  by  using  GPS  observations  (see  Fig.  14).  The  origin  is  the  center  of  the  Earth,  and  the  car  is  originally 
located  at  the  surface  of  the  Earth — (0,  Re)  in  an  Earth  centered  coordinate,  where  Re  =  6357  km  is  the  radius 
of  the  Earth — which  is  assumed  to  be  unknown.  The  velocity  can  be  read  from  the  velocimeter,  and  has  a 
constant  true  value  of  70  mph,  or  31.3  m/s.  The  velocity  is  assumed  to  have  a  Gaussian  measurement  error  with 
covariance  matrix  Ev  =  a l  diag[  1  1],  where  av  is  regarded  as  an  adjustable  parameter.  The  direction  of  the 
velocity  changes  in  time.  The  car  is  equipped  with  a  GPS  receiver.  Four  visible  satellites  (ns  =  4)  are  located 
at  the  altitude  of  20200  km  above  the  Earth  surface  with  initial  angles  7r/5,  77t/24,  47t/7,  27t/3  (in  the  Earth 
centered  coordinate).  The  satellites  orbit  the  Earth  at  a  period  of  12  hours.  The  receiver  on  the  car  can  receive 
GPS  signal  from  each  satellite  at  the  frequency  of  1  Hz  (A t  =  tk  -  tk- 1  =  1  s),  from  which  the  distances  from 
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Figure  14:  Setup  of  the  numerical  problem,  (a)  The  tracks  of  the  satellites  and  the  car  for  the  simulation  time 
(500  seconds).  The  satellites  move  counterclockwise,  (b)  The  track  of  the  car,  where  the  origin  is  shifted  to  the 
original  position  of  the  car. 


Figure  1 5:  For  the  two-dimensional  positioning  problem,  dependence  of  the  standard  deviations  of  the  samples  ox 
for  the  two  coordinates  on  the  standard  deviation  of  the  velocity  ov.  The  standard  deviation  of  the  pseudoranges 
is  op  —  2.5,  the  number  of  samples  is  N  =  1000,  and  h  =  0.3.  Symbols  are  obtained  from  numeric  simulations, 
where  each  data  is  the  average  of  10  runs,  and  100  different  time  steps  for  each  run  are  used.  The  curves  are  from 
the  our  theoretical  prediction  Eq.  (37). 
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Figure  16:  A  comparison  of  the  errors  of  the  first  coordinate  x\  of  the  position  prediction  by  velocity  information 
only  (a),  by  least  squares  of  GPS  distances  with  outliers  only  (b),  and  by  our  strategy  (c).  The  parameter  values 
used  in  the  simulation  are  av  —  0.5,  ap  =  2.5,  cq  =  0.7,  c\  =  4.2,  C2  =  7,  N  =  1000,  and  h  =  0.3. 


the  satellites  to  the  car  Pk  (pseudoranges)  can  be  measured.  Assuming  that  the  receiver  has  no  clock  offset  and 
the  satellites  are  not  correlated,  we  can  write  the  covariance  matrix  of  the  pseudoranges  as  Ep  =  a 2  7n,,  where 
Ina  is  the  identity  matrix  of  order  ns.  A  motion  model  which  is  linear  in  the  state  dynamics  and  nonlinear  in  the 
measurements  is  [59] 


x(t)  =  x{tk- 1)  +  v(tk- 1)  •  A t, 
y(t)  =  9[x{t)]  +  cp(0, 

where  y  is  the  pseudorange  measurements  y  =  [P1  P 2  •••  PUa]T .  The  measurement  function  is  g(x)  = 
[Rl  R2  •  •  •  Rna]T,  where  Ri  =  ||XJ  —  x||  is  the  Euclidean  distance  from  the  car’s  position  x  to  the  j' th 
satellite  A7,  and  ep  is  the  pseudorange  observational  noise.  The  covariance  matrices  for  v  and  ep  are  Et,  and  Ep 
respectively. 

Figure  15  shows  the  dependence  of  Ex  on  Ev  when  Es  is  given,  which  can  be  obtained  numerically  from  the 
distribution  of  the  pseudoranges  (Ep).  The  symbols  are  obtained  from  direct  simulations,  the  curves  are  from  our 
theory  Eq.  (37).  They  agree  quite  well. 

To  test  the  robustness  of  our  Monte-Carlo  strategy,  we  added  15  random  outliers  of  30  meters  to  the  measured 
pseudoranges  of  the  second  satellite  in  a  time  span  of  500  seconds  with  measurement  frequency  1  Hz.  The  result 
of  position  estimation  is  shown  in  Fig.  16.  Three  cases  are  presented  for  comparison.  Figure  16(a)  shows 
the  prediction  error  of  x\  where  the  initial  position  is  known  and  the  measured  velocity  has  standard  deviation 
ov  —  0.5.  There  is  a  systematic  drift  of  errors.  Figure  16(b)  shows  the  prediction  error  of  the  least-squares  method 
[72]  if  only  the  measured  GPS  pseudoranges  are  available  (standard  deviation  op  =  2.5  for  each  satellite).  Figure 
16(c)  is  the  estimation  error  from  our  proposed  Monte-Carlo  strategy,  which  apparently  contains  no  systematic 
error  and  exhibits  much  smaller  statistical  errors. 

Next,  we  studied  the  cases  with  non-Gaussian  noise  in  the  GPS  pseudorange  observations.  This  might  be 
the  case  when  the  car  is  moving  in  a  city  or  in  the  forests,  where  GPS  signal  may  be  blocked  by  the  buildings 
or  the  trees  and  causes  difficulty  to  distinguish  multipath  signal  from  the  original  signal,  which  may  introduce 
systematic  biases  [2].  Furthermore,  because  of  the  complexity  of  the  environment,  at  certain  moments  the  original 
signal  may  be  unavailable.  We  assume  that  the  distribution  of  noise  in  the  pseudoranges  consists  of  two  Gaussians 
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Figure  17:  A  comparison  of  the  errors  of  the  second  coordinate  X2  of  the  position  prediction  only  by  least  square 
of  GPS  distances  with  outliers  (a),  by  regularized  sequential  Monte-Carlo  simulation  without  the  robust  strategy 
(b),  and  by  our  strategy  (c).  The  parameter  values  used  in  the  simulation  are  av  =  1,  ap  =  2.5,  cq  =  1.5, 
ci  =  4.2,  C2  =  8,  N  =  1000,  and  h  =  0.3. 


with  different  mean  values.  The  probability  density  function  is: 

1  _  X2  |  _  (J  — Jp)2 

/(*)  =  b  /7T~  e  *%  +  (!-  b)-7^=—e  20 2  -  (39) 

V27T(Ti  V27TCT2 

where  6  is  a  weight  factor  and  we  take  b  =  0.6  in  our  simulation.  Other  parameters  are  o\  =  ov  —  2.5,  ^2  =  1, 
xo  =  8. 

Again,  to  test  the  robustness  of  our  algorithm,  we  added  outliers  to  the  GPS  pseudorange  observations:  20 
outliers  of  40  meters  are  added  to  satellite  2  randomly.  Figure  17  compares  the  errors  in  position  estimation 
by  three  methods:  the  least-squares  estimation  from  the  GPS  pseudoranges  with  outliers  (a),  the  estimation  by 
the  regularized  sequential  Monte-Carlo  scheme  (b),  and  by  our  robust  sequential  Monte-Carlo  scheme  (c).  The 
least-squares  estimation  from  the  pseudoranges  has  large  errors  and  can  have  systematic  deviations  (the  average 
of  the  error  is  not  zero),  as  shown  in  Fig.  17(a).  The  regularized  sequential  Monte-Carlo  scheme  can  remove 
these  systematic  deviations  caused  by  the  non-Gaussian  distribution  but  is  affected  by  the  outliers,  which  can 
be  seen  from  the  spikes  in  Fig.  17(b).  Our  robust  sequential  Monte-Carlo  scheme  can  recover  from  both  the 
systematic  deviations  and  the  outliers  [Fig.  17(c)].  In  fact,  the  average  absolute  value  of  the  errors  can  be  30% 
smaller. 

The  current  robust  scheme  deals  with  observational  outliers.  If  there  are  dynamical  outliers — e.g.  the  outliers 
appear  in  v — the  current  scheme  needs  to  be  modified  to  cope  with  the  problem.  Observations  after  such  an  event 
will  be  needed  to  identify  an  outlier. 

3.3.6  Discussions 

In  conclusion,  we  have  obtained  a  self-consistent  equation  for  the  estimation  precision  of  the  Bayesian  inference 
in  terms  of  the  dynamical  noise  and  the  observational  noise  levels.  The  equation  may  provide  insights  into 
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designing  improved  sequential  Monte-Carlo  simulations  with  higher  precision.  We  also  proposed  a  strategy  to 
deal  with  sudden,  large  disturbances  that  are  inevitable  in  physical  observations.  The  effectiveness  of  our  method 
has  been  tested  numerically.  While  we  used  a  standard  kinematic  GPS  positioning  problem  to  demonstrate  the 
working  of  our  approach,  it  is  general  and  applicable  to  signal-processing  problems  where  outliers  are  present 
and  are  to  be  mitigated.  Sequential  Monte-Carlo  simulations  have  begun  to  be  used  widely  in  various  estimation 
problems  in  science  and  engineering.  Our  contribution  provides  a  robust  strategy  for  improving  the  estimation 
precision  when  experimental  observations  are  nonstationary  or  even  temporally  interrupted. 
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6  Interactions/Transitions 

There  was  an  active  collaboration  with  Drs.  S.  Mohamod  and  J.  Murchison  at  Eglin  AFB.  The  activities  resulted 
in  one  joint  paper  published  in  Journal  of  Geodesy. 

During  the  project  period,  the  PI  gave  two  invited  talks  on  GPS,  in  addition  to  talks  at  annual  AFOSR  Test 
and  Evaluation  Workshops. 

1 .  “Performance  of  integer  least-squares  method  for  ambiguity  resolution  with  GPS  signals  in  the  presence  of 
noise,”  Seminar,  Tamasek  Laboratory  and  Department  of  Physics,  National  University  of  Singapore,  July 
20,  2004. 
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